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ABSTRACT 

We report the discovery of HAT-P-31b, a transiting exoplanet orbiting the V=11.660 dwarf star 
GSC 2099-00908. HAT-P-31b is the first HAT planet discovered without any follow-up photometry, 
demonstrating the feasibility of a new mode of operation for the HATNet project. The 2.17 Mj, 
1.1 Rj planet has a period Pb = 5.0054 days and maintains an unusually high eccentricity of ei, = 
0.2450 ± 0.0045, determined through Keck, FIES and Subaru high precision radial velocities. Detailed 
modeling of the radial velocities indicates an additional quadratic residual trend in the data detected 
to very high confidence. We interpret this trend as a long-period outer companion, HAT-P-31c, of 
minimum mass 3.4 Mj and period > 2.8 years. Since current RVs span less than half an orbital period, 
we are unable to determine the properties of HAT-P-31c to high confidence. However, dynamical 
simulations of two possible configurations show that orbital stability is to be expected. Further, if 
HAT-P-31c has non-zero eccentricity, our simulations show that the eccentricity of HAT-P-31b is 
actively driven by the presence of c, making HAT-P-31 a potentially intriguing dynamical laboratory. 

Subject headings: planetary systems — stars: individual (HAT-P-31, GSC 2099-00908) techniques: 
spectroscopic, photometric 



1. INTRODUCTION 

Transiting extrasolar planets provide invaluable in- 
sight into the nature of planetary systems. The op- 
portunities for follow-up include spectroscopic infer - 
ence of an exoplanet's atmosphere dTinetti et al. [2007), 
searches for dy namical variations ( Agol et al. 120051 : 
IKipping "| [2009a, b), and characterizing the orbital ele- 
ments ( Winn et al.ll20lil) . Multi-planet systems in par- 
ticular offer rich dynamical interactions and their fre- 
quency is key to understanding planet formation. 

The Hu ngarian- made Aut omated Telescope Network 
(HATNet; Ba kos et al.ll200l survey for transiting exo- 
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planets (TEPs) around bright stars operates six wide- 
field instruments: four at the Fred Lawrence Whipple 
Observatory (FLWO) in Arizona, and two on the roof 
of the hangar servicing the Smithsonian Astrophysical 
Observatory's Submillimeter Array, in Hawaii. Since 
2006 . HATNet has anno unced and published 30 TEPs 
(e.g. I Johnson et al.l |20TTT ). In this work, we report our 
thirty-first discovery, around the relatively bright star 
GSC 2099-00908. In addition, a long-period companion 
is detected through detailed modeling of the radial ve- 
locities, although no transits of this object have been 
detected or are necessarily expected. 

In Section [2] we summarize the detection of the pho- 
tometric transit signal and the subsequent spectroscopic 
observations of HAT-P-31 to confirm the planet. In 
Section [3j we analyze the data to determine the stellar 
and planetary parameters. Our findings are discussed in 
Section |U 

2. OBSERVATIONS 

As described in detail in several pr evious papers (e.g. 
iBakos et al.ll2010HLatham et alj|2009fl . HATNet employs 
the following method to discover transiting planets: 1. 
Identification of candidate transiting planets based on 
HATNet photometric observations. 2. High-resolution, 
low-S/N "reconnaissance" spectra to efficiently reject 
many false positives. 3. Higher-precision photometric 
observations during transit to refine transit parameters 
and obtain the light curve derived stellar density. 4. 
High- resolution, high-S/N "confirmation" spectroscopy 
to detect the orbital motion of the star due to the planet, 
characterize the host star, and rule-out blend scenarios. 

In this work, step three is omitted and this is the 
biggest difference to the usual HATNet analysis. The 
detection, and thus verification, of an exoplanet can be 
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made using step four alone. Indeed, the majority of ex- 
oplanets have been found in this way. Step 1 clearly al- 
lows us to intelligently select the most favorable targets 
for this resource-intensive activity though. Step 2 follows 
the same logic. Step 3 is predominantly for the purposes 
of characterizing the system, and therefore its omission 
does not impinge on the planet detection. In some cases, 
follow-up photometry is used to confirm marginal HAT- 
Net candidate detections as well, but this is not the case 
for the discovery presented in this work. An additional 
check on step 4 is that the derived ephemeris is consistent 
with that determined photometrically. 

We did not obtain high precision photometry for this 
target as the transits were not observable from our usual 
site of choice, FLWO, until at least May 2012. This 
is because the transiting planet has a near-integer pe- 
riod and the time of transit minimum has now phased 
into the day-light hours in Arizona (i.e. unobservable) . 
Rather than wait until this time, we have decided to re- 
lease this confirmed planet detection to the community 
so that follow-up photometry may be conducted at other 
sites. 

The principal consequence of not having any follow-up 
photometry is that the obtainable precision of the transit 
parameters is reduced. In the case of HAT-P-31b, this 
means that the light curve derived density was less pre- 
cise than that determined spectroscopically. In practice 
then, we reverse the usual logic and instead of apply- 
ing a prior on the stellar density from the light curve, 
we apply a prior on the light curve from the stellar den- 
sity. One can see that the decision on this will vary from 
case to case depending upon the transit depth and target 
brightness. 

Another issue is that in the past HAT analyses have 
used the HATNet photometry for measuring the plane- 
tary ephemeris, P and r, and little else. All other param- 
eters could be more precisely determined from the follow- 
up photometry. As a consequence, w e used the Externa l 
Parameter Decorrelation (EPD; see iBakos et al.l l2010f) 
and T rend Filtering Algorithm (TFA: see IKovacs et al.1 
2005) techniques to correct the HATNet photometry, 
which are known to attenuate the apparent transit depth 
by a small amount. This was usually accounted for by 
including an instrumental blending factor, -Bi ns t, in the 
HATNet data, which could be determined by comparing 
the ratio of the HATNet apparent depth and the follow- 
up photometry depth. Without follow-up photometry, 
-Bi ns t would be unconstrained and so an estimation of the 
planetary radius would be impossible. To avoid this, we 
employ the more computationally demandi ng and sophis- 
ticate d technique of reconstructive TFA ([Kovacs et al.l 
120051) . Reconstructive TFA does not attenuate the tran- 
sit depth significantly and thus offers a way of avoiding 
a free -Bi ns t term. Our previous experience with the two 
modes of TFA sup port this. For exam ple, HAT-P-15b's 
TFA photometry (|Kovacs et al.l I201QT) was found to re- 
quire a blending factor of B mst = 0.71 ± 0.07. In con- 
trast, the reconstructive TFA for the same planet causes 
-Bi ns t = 0.95 ± 0.04. We find no instance in any previous 
implementation of reconstructive TFA where -Binst would 
depart from unity by more than 2-cr. We will therefore 
use the reconstructive TFA in this work and conserva- 
tively double all uncertainties relating to the depth and 
radius of the planet. 



In this paper, we will consequently show that HATNet 
photometry alone is sufficient to constrain the system 
properties and that future work may not always require 
step 3 (i.e. follow-up photometry). 

In the following subsections we highlight specific de- 
tails of this procedure that are pertinent to the discovery 
of HAT-P-31b. 

2.1. Photometric Detection 

The transits of HAT-P-31b were photometrically de- 
tected with a combined confidence of 6.2-ct using the 
HAT-5 telescope in Arizona and the HAT-8 telescope 
in Hawaii. The region around GSC 2099-00908, a field 
internally labeled as 241, was observed between 2007 
March and 2007 July, whenever weather conditions per- 
mitted. In total, we gathered 9205 exposures of 5 min- 
utes at a 5.5 minute cadence. 769 of these images were 
rejected by our reduction pipeline because they produced 
bad photometry for a significant fraction of stars. A typ- 
ical image is found to contain approximately 48,000 stars 
down to Ic ~ 14. For the brightest stars in the field, the 
photometric precision per-image was 3 mmag. 




0.8 1 1.2 
Frequency [d ] 



Fig. 1.— Box-fitting Least Squares (BLS; [Ko7acs et al.1120013) 
periodogram of HATNet photometry. The arrow marks the true 
transit signal and the other peaks are inter pretted to be a liases . 
The y-axis denotes Signal Residue (SR), sec Kovacs ct al. (2002) 
for details of the definition. 

Standard photometric procedures were used to cali- 
brate the HATNet frames and then these calibrated im- 
ages were subj ected to star detec tion and astrometry, 
as described in lPal fc B akos (2006). Aperture photome- 
try was performed on each image at the stellar centroids 
derived from the Tw o Micron All Sky Survey (2MASS; 
ISkrutskie et al.ll2006D catalog and the individual astro- 
metric solutions. The resulting light curves were decor- 
related (cleaned of trends) using the Externa l Parame- 
ter Decorrelation (EPD; see lBakos et al.ll2010D technique 
in "const ant" mode and the Trend Filtering Algorithm 
(TFA; see IKovacs e1~aTl l2005h . 

The light curves were searched for transits us- 
ing t he Box-fitting Least-Squares (BLS; IKovacs et al.l 
2002) method. We detected a significant signal in 
the light curve of GSC 2099-00908 (also known as 
2MASS 18060904+2625359; a = 18 h 06 m 09.00 s S = 
+26°25'36.0"; J2000; V=11.660 IDroeee et all I2006T) . 
with an apparent depth of ~ 5.1 mmag, and a period 
of P = 5.0050 days. The BLS periodogram is shown 
in Figure [TJ The drop in brightness had a first-to- 
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Fig. 2.— Unbinned light curve of HAT-P-31 including all 8436 
instrumental Iq band 5.5 minute cadence measurements obtained 
with the HAT-5 and HAT-8 telescopes of HATNet (see the text for 
details), and folded with the period P = 5.005425 days resulting 
from the global fit described in Section [3} . T he so lid line shows 
our transit model fit to the light curve (Section 13, 2|l . The bottom 
panel shows a zoomed-in view of the transit, the filled circles show 
the light curve binned in phase with a bin-size of 50 points. 



TABLE 1 

HATNet differential photometry of HAT-P-31 



BJD 


Mag (EPD) a 


Mag (TFA) b 


°"Mag 


(2,400,000+) 








54178.0007000 


11.4164 


11.4203 


0.0075 


54178.0045362 


11.4267 


11.4273 


0.0073 


54178.0083829 


11.4250 


11.4317 


0.0075 


54178.0122220 


11.4177 


11.4275 


0.0069 


54178.0160619 


11.4244 


11.4220 


0.0072 



Note. — This table is available in a machine-readable form 
in the online journal. A portion is shown here for guidance 
regarding its form and content. 

a These magnitudes have been subjected to the EPD proce- 
dure. 

b These magnitudes have been subjected to the EPD and TFA 
procedures. 

last-contact duration, relative to the total period, of 
q = 0.0432 ± 0.0038, corresponding to a total duration 
of 5.187 ± 0.462 hr. Due to the lack of follow-up pho- 
tometry, the HATNet photometry was re-processed with 
more computationally expensive reconstructive TFA, as 
discussed in $21 an d shown in Figure [5] The EPD and 
reconstructive TFA corrected photometry is provided in 
Table [U 

2.2. Reconnaissance Spectroscopy 



High-resolution, low-S/N reconnaissance spectra were 
obtained for HAT-P-31 using the Tillinghast Reflec- 
tor Echelle Spectrograph (TRES; iFuresz et all 120081) 
on the 1.5 m Tillinghast Reflector at FLWO, and the 
echelle spectrograph on the Australian National Univer- 
sity (ANU) 2.3 m telescope at Siding Spring Observatory 
(SSO) in Australia. The two TRES spectra of HAT-P- 
31 were obtained, reduced and analyzed to measure the 
stellar effective temperature, surface gravity, projected 
rotation velocity, and RV via cross-correlation against a 
library of synthetic template spectra. The r eduction and 
analys is pr ocedure has been des cribed by IQuinn et al.l 
(2010) and lBuchhave et al.l (|2010f >. A total of 14 spectra 
of HAT-P-31 were obtained with the ANU 2.3 m tele- 
scope. These data were collected, reduced and analyzed 
to measure the RV via cross-correlation against the spec- 
trum of a RV standar d star HP 22331 1 following the 
procedure described bv lBekv et al.l (|2011|) . The resulting 
measurements from TRES and the ANU 2.3 m telescope 
are given in Table [2] 

These observations revealed no detectable RV varia- 
tion at the lkms -1 precision of the observations. Addi- 
tionally the spectra are consistent with a single, slowly- 
rotating, dwarf star. 

2.3. High Resolution, High S/N Spectroscopy 

We proceeded with the follow-up of this candidate by 
obtaining high-resolution, high-S/N spectra to charac- 
terize the RV variations, and to refine the determina- 
tion of the stellar parameters. For this we used the 
HIR ES instrument ifVogt et al.l [l99l with the iodine- 
cell (jMarcv fc Butler! I1992D on the Kec k I telescope 
the H igh-Dispersion Spect rograph (HPS: | Noguchi et al.l 
12002ft with the iodine-cell (|Sato et alJl2002ft on the Sub- 

aru telescope, and the FIbr-fed Echelle Spectrograph 
(FLES) on the 2.5 m Nordic Optical Telescope (NOT; 
iPiupvik fc AnderserJ l2010). Tabled summarizes the ob- 
servations. The table also provides references for the 
methods used to reduce the data to relative RVs in the 
Solar System barycentric frame. The resulting RV mea- 
surements and their uncertainties are listed in Table |H 
The different instrumental uncertainties arise from dif- 
ferent slit widths, exposure times and seeing conditions. 
The period-folded data, along with our best fit described 
below in Section [3j are displayed in Figure [3] 

One false-alarm possibility is that the observed radial 
velocities are not induced by a planetary companion, 
but are instead caused by distortions in the spectral line 
profiles due to contamination from a nearby unresolved 
eclipsing binary. This hypothesis may be interrogated 
by examining the spectral line profiles for c ontamination 
from a nearby unresolve d eclipsing binary (jQueloz et al.1 
I200H iTorres et al.l 120071 ). A bisector analysis based on 
the Keck spectra was performed as described in §5 of 
iBakos et al.l ([2007) . The resulting bisector spans, plot- 
ted in Figure [3[ show no significant variation, and are 
not correlated with the RVs, indicating that this is a real 
TEP system. 

In the same figure, one can a lso see the S index 
(|Vaughan. Preston fc Wilson! 119781) . which is a quanti- 
tative measure of the chromospheric activity of the star 
derive d from the flux in the co res of the Ca II H and K 
lines (jlsaacson fc Fischer] 12010ft . Following iNoves et al.l 
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TABLE 2 

Summary of reconnaissance spectroscopy observations of HAT-P-31 



Instrument 


Date 




Number of 




l°g(ff* [cgs]) 


v sini 


7RV a 








Spectra 


[K] 




[kms -1 ] 


[kms -1 ] 


TRES 


2009 Jul 


05 


1 


6000 


4.0 


4 


-2.342 


TRES 


2009 Jul 


07 


1 


6000 


4.0 


4 


-2.300 


ANU 2.3 m 


2009 Jul 


14 


5 








-8.07 ±0.35 


ANU 2.3 m 


2009 Jul 


17 


5 








-7.13 ±0.40 


ANU 2.3 m 


2009 Jul 


18 


2 








-7.64 ±0.44 


ANU 2.3 m 


2009 Jul 


19 


2 








-7.54 ±0.49 


a The mean heliocentric RV 


of the target. 


Systematic differences between the velocities from 



the two instruments arc consistent with the velocity zero-point uncertainties. For the ANU 2.3 m 
observations we give the weighted mean of the observations and the uncertainty on the mean 
for each night. Note that the systematic difference of — 5.3kms _1 between the ANU 2.3m and 
TRES observations is similar to the difference of —5.1 kms -1 found between these same two 
instruments bv lBeky eFaXI d2011fl for HAT-P-27. 



TABLE 3 

Summary of high resolution, high S/N spectroscopy 

OBSERVATIONS OF HAT-P-31 



Instrument 


Date 


Number of 


Reduction 




Range 


RV obs. 


Reference 


HDS 


2009 Aug 08 - 2010 May 24 


25 


1 


HIRES 


2010 Feb 24 - 2010 Jul 3 


9 


2 


FIES 


2009 Oct 6 - 2009 Oct 11 


6 


3 



References. — 1: ISato et all 1 12009) . 2: IButler et al] U996I 1. 3: 
IBuchhave et al.l II2010T) 

a The mean heliocentric RV of the target. Systematic differences between 
the velocities from the two instruments arc consistent with the velocity zero- 
point uncertainties. 



we find that HAT-P-31 has an activity index 
logi?^ K = —5.312, implying that this is a very inactive 



star. 



3. ANALYSIS 



The analysis of the HAT-P-31 system, including deter- 
minations of the properties of the host star and planet, 
was carried out i n a similar fashion to previous HATNet 
discoveries (e.g. iBakos et al.l l2010h . Below, we briefly 
summarize the procedure and the results for the HAT- 
P-31b system. 

3.1. Properties of the Parent Star 

Stellar atmospheric parameters were measured from 
our template Keck/HIRES spectrum using the Spec- 
troscopy Made Easy (SME; iValenti fc Piskunovl 11996ft 
analysis package, and the atomic line database of 
IValenti fe Fischer! (T2005h ■ SME yielded the following 
values and uncertainties: effective temperature T e ff* = 
6065±100K, metallicity [Fe/H] = +0.15 ± 0.08dex, and 
stellar surface gravity log = 4.261q;i3 ( c S s )' projected 
rotational velocity usini = 0.5 ± 0.6 kms -1 . 

The above atmospheric p arameters are t hen combined 
with the Yonsei-Yale (YY) (|Yi et alJl20"ol series of stel- 
lar evolution models to determine other parameters such 
as the stellar mass, radius and age. The results are listed 
in Table We find that the star has a mass and radius 
of M* = 1.218±g;g||M and i?* = 1.36^;^ i? Q , and 

an estimated age of 3.17^5'jj 1 Gyr. 

For previous HATNet planets (e.g. IBakos et al.ll2010D 
we used the normalized semimajor axis, a/i?* which is 



closely related to p*, the mean stellar density, and is 
determined from the analysis of the light curves and RV 
curves, to obtain an improved estimate of logg*, which 
is then held fixed in a second SME iteration. In this 
case, because we only have the relatively low-precision 
HATNet light curve, a/i?* is poorly constrained, and we 
instead opt to use the SME determination of log g+ rather 
than a/Ri, as a luminosity indicator. Figure 3] shows the 
location of the star in a diagram of logg* versus T e ff*, 
together with the model isochrones. For comparison we 
also show the relatively poor constraint on log 5* that is 
imposed by a/R*. 

As an additional check on the stellar evolution mod- 
eling, we note that HAT-P-31 has a measured near- 
infrared color of J — K = 0.364 ± 0.034, which we 
have taken , from 2MASS (|Skrutskie et al.ll2006| ) using the 
iCarpenterl ()2001[ ) transformation to the ESO photomet- 
ric system. This is within 2-cr of the predicted value 
from the isochrones of J — K = 0.31 ± 0.11. The dis- 
tance listed in Table [3] is calculated by comparing the 
observed K magnitude (taken from 2MASS and trans- 
formed to ESO) to the absolute K magnitude from the 
stellar models. 

3.2. Global Modeling of the Data 
3.2.1. Photometry 

In previous HATNet papers, we have used a simplified 
model for the transit light curve of the HATNet data. For 
HAT-P-31, no precise photometry exists and thus we fit 
the HATNet da ta using a more s ophis ticated quadratic 
limb darkening iMandel fe Agoll (|2002l ) algorithm with 
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TABLE 4 

Relative radial velocities, bisector spans, and activity index 
measurements of hat-p-31. 



BJD a 


RV b 


°-rv c 


BS 


CTBS 


S d 


Instrument 


[2,454,000+] 


[rns^ 1 ] 


[ms" 1 ] 


[ms- 1 ] 


[ms ] 






1051.87826 


—30.72 


6.99 








Subaru 


1 n t: i oon a c 

1051.8894b 


—27.97 


6.90 








Subaru 


1 n tr i r\r\r\t~"~7 

1U51.9UUD7 


—44.27 


6.56 








Subaru 


iriKO Q1 /ITI 


— 163.76 


9.25 








Subaru 


1 nrn no 000 

1052.93238 


oon on 

— 200.89 


6.63 








Subaru 


lU5z.947Uo 


—202.52 


7.29 








Subaru 


IflCQ OH /I Q Q 

1Uoo.oU4oo 


— 159.83 


7.09 








Subaru 


105a. 81555 


— 166.62 


7.47 








Subaru 


1053.8267b 


— 163.36 


7.95 








Subaru 


1 n rz on rz a a 

1053.89544 


— 145.48 


7.08 








Subaru 


1 n et nr\cc tz 

1053.90bb5 


— 140.21 


6.46 








Subaru 


1053.91786 


— 136.05 


6.39 








Subaru 


1111 OOylOO 


71.48 


9.40 








FIES 


1 1 1 O QQ "7Tr\ 

1112.66 J iv 


— 115.00 


10.70 










1110 r 

lllo.o4Ulo 


—225.79 


10.50 








FIES 


111/1 O /I wo 

1114. 34798 


12.92 


10.70 








FIES 


1115.38b54 


234.94 


10.80 








FIES 


iiii^ oo/im 

lllo.oo4Ul 


101.33 


24.70 








FIES 


1 nrn 1 n A C tZ 

12oz.lU4oo 






—2.81 


1.78 


0.1280 


Keck 


lzoz. 11375 


—34.70 


2.65 


5.07 


1.72 


0.1190 


Keck 


lz8o.U98/ f 


153.49 


2.56 


2.47 


1.48 


0.1270 


rveck 


1QQO OOC7Q 

1660. yuo (6 


— 199.81 


11.70 








Subaru 


1 000 m 1 nn 

1338.911UU 


— 199.43 


11.83 








Subaru 


l0OO.9l0Z0 


— 196.43 


10.73 








Subaru 


1338. 91953 


— 188.53 


11.15 








Subaru 


1 o on ocrn-icr 

1339.85945 


129.95 


16.85 








Subaru 


1 oin f7or\A 

1339.87200 


140.14 


11.90 








Subaru 


1339.88b69 


151.20 


9.63 








Subaru 


1339.95891 


180.56 


7.45 








Subaru 


1339.97360 


179.56 


7.50 








Subaru 


1339.98489 


186.81 


8.74 








Subaru 


1341.08414 


167.01 


13.93 








Subaru 


1341.09189 


158.61 


13.28 








Subaru 


1341.10303 


151.27 


9.77 








Subaru 


1343.00801 


-167.85 


2.30 


10.12 


1.17 


0.1350 


Keck 


1372.86477 


-139.73 


2.11 


-3.69 


1.30 


0.1230 


Keck 


1374.99015 


177.54 


1.84 


0.74 


1.16 


0.1240 


Keck 


1375.86758 


207.17 


1.90 


-3.49 


1.11 


0.1240 


Keck 


1378.79885 


-228.26 


1.90 


-4.23 


1.29 


0.1230 


Keck 


1380.79961 


212.57 


1.93 


-4.20 


1.12 


0.1220 


Keck 



Note. — For the iodine-free template exposures there is no RV measurement, 
but the BS and S index can still be determined. 

a Barycentric Julian dates throughout the paper arc calculated from Coordinated 
Universal Time (UTC). 

k The zero-point of these velocities is arbitrary. An overall offset 7 re i fitted to these 
velocities in Scction l3.2l has not been subtracted. 

c Internal errors excluding the component of astrophysical/instrumcntal jitter con- 
sidered in Scction l3.2l 

^ Relative chromospheric activity index, not calibrated to the scale of 
IVaughan. Preston fe Wilson! 119781) . 

li mb dark e ning c oefficients interpolated from the tables 
by [Claret (2004). One caveat is that the instrumental 
blending factor, £?i ns t, is unknown as discussed earlier in 
§21 We point out that experience with previous HAT- 
Net planets suggests Bi nst is within 2-cr of unity for light 
curves processed using reconstructive TFA in all cases 
and thus can be accounted for by conservatively dou- 
bling the uncertainties on p and Rp. Further support 
for a B mst factor not greatly different from unity come 
from the fact HAT-P-31 is fairly isolated and there are 
no neighbors in 2MASS or a DSS image which contribute 
significant flux to the HATNet aperture. 

Due to the low-precision photometry, the stellar den- 
sity canno t be determined to high precisio n using the 
method of S eager fc Mallen-Ornelasl ()2003l ) and in fact 
spectroscopic estimates were found to be more precise. 



However, we can reverse this well-known trick by imple- 
menting a Bayesian prior in our fitting process for the 
stellar density. 

We use the spectroscopically determined stellar den- 
sity from §3.11 as a prior in our fits. Since the period 
of the transiting planet is well constrained for even low 
signal-to-noise photometry, reasonably precise estimates 
for P and p* are possible. With these two constrained, 
(db/R*) is therefore also constrained. The transit light 
curve is essentially characterized by four parameters, t&, 
8b, bb and (db/R*) and thus one of these parameters is 
constrained by the combination of the transit times and 
the stellar density prior alone. 

The photometry which is fitted in the global model- 
ing is corrected for instrumental systematics through the 
EPD and reconstructive TFA correction procedures prior 
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Fig. 3. — first Row; Keck/HIRES (squares), Subaru (circles) and FIES (triangles) RV measurements for HAT-P-31, along with our 
best-fit 2-planet model (see Table [Jj. The center-of-mass velocity has been subtracted. Second Row: Same as top panel except the RV 
model of the inner planet has been subtracted from the data and the model, revealing the orbit of the outer planet. The x 2 of the best-fit 
2-planet model is 34.8 for 39 data points (rms of 9.60 m/s) indicating a stellar jitter at or below the measurements errors. Third Row: 
Residuals from our best-fit model. Lower left: RV measurements phased to the orbital periods of the inner planet (left). The quadratic 
trend has been removed. Lower right: Upper panel shows the bisector spans (BS), with the mean value subtracted, phased at the period 
of the inner planet. Lower panel shows relative chromospheric activity index S measured from the Keck spectra, phased at the period of 
the inner planet. Note the different vertical scales of the panels. Observations shown twice are represented with open symbols. 
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Fig. 4. — Model isochrones from Yi ct al. (2001) for the measured 
metallicity of HAT-P-31, [Fe/H]= 0.15, and ages from 1 to 14Gyr 
in steps of 1 Gyr (gray-scale lines, left to right). The adopted 
values of T e g t and log determined from the SME analysis, are 
shown together with their 1-a and 2-er confidence ellipsoids (filled 
circle, and bold solid lines). The less-tight 1-a and 2-<r constraints 
imposed by a/R+ are also shown (open triangle, and bold dotted 
lines). 



eccentric transiting planet with a quadratic trend in the 
RVs was the preferred model. The pivot point (i p i vo t) of 
the polynomial models, including the drift and quadratic 
trends, was selected to be the weighted mean of the radial 
velocity time stamps. 

The most likely physical explanation for a quadratic 
trend is a third body in the system, described by a Ke- 
plerian model. Indeed, the Keplerian model provides an 
improved % 2 for a circular orbit and then again for an 
eccentric orbit but the extra degrees of freedom penalize 
our model selection criterion. We also found that these 
models were highly unconstrained and convergence in the 
associated fits was unsatisfactory. An illustration of the 
lack of convergence is shown in Figure [6l Therefore, we 
will adopt the quadratic model in our final reported pa- 
rameters in Table [7J 

The quadratic model may be used to infer some physi- 
cal parameters of the third planet. To make some mean- 
ingful progress, we will assume the outer planet is on 
a circular orbit. One may compare the quadratic and 
Keplerian model descriptions via: 



TABLE 5 
Stellar parameters for HAT-P-31 



Parameter 



Value 



Source 



Spectroscopic propcrtii 



T cfu IK] 

[Fe/H] 


6065 ± 100 


SME a 


0.15 ±0.08 


SME 


v sin i [km s 1 ] . . . . 


0.5 ±0.6 


SME 


"mac [kms~ 1 ] b . . . 


4.47 


SME 


Vjvic [kms _1 ] b 


0.85 


SME 


7RV [kms- 1 ] 


-2.40 ±0.03 


TRES 



Photometric properties 

V [mag] 11.660 

V-I c [mag] 0.67 ±0.17 

J [mag] 10.423 ± 0.023 

H [mag] 10.128 ±0.027 

Ks [mag] 10.083 ± 0.021 

Derived properties 



TASS 

TASS 

2MASS 

2MASS 

2MASS 



M* [Mq] 


1 ois+U' ™ 
1 - zlo -0 .063 


YY+SME c 


R* [Rq] 


^"-o.is 


YY+SME 


l°g(9* [cgs]) 


4 26+° 11 


YY+SME 


L* [£©1 


9 9 ,+1.01 
^• ZO -0.58 


YY+SME 


My [mag] 


, qi+0.34 


YY+SME 


M K [mag,ESO] . . . 


2.64 ± 0.30 


YY+SME 


Age [Gyr] 


, I 7 +0.70 


YY+SME 


Distance [pc] 




YY+SME 



a SME — "Spectroscopy Made Easy" package for the anal- 
ysis of high-resolution spectra (Valcnti & Piskunov 1996). 
k Assumed quantity based upon derived spectral type. 
c YY+SME = Based on the YY isochrones lIYi et al.H200in 
and the SME results. 



to performing the fit (see 
for details). 



i CO] and IBakos et all lIMoh 



3.2.2. Radial velocities 

For the radial velocity fits, we found a single planet fit 
gave a very poor fit to the observations (x 2 = 204.2 for 
40 RV points) and quickly appreciated some kind of sec- 
ond signal was present. Exploring different models, such 
as Trojan offsets, polynomial time trends and outer com- 
panions, (see Table |6] for comparison), we found that an 



RV^ uad = 7' + 7(t - ipivot) + 0.5 7 (i - Vvot) 2 (!) 
'27r(t-r c )\ 



RV^ ep = 7 '-^cSin 



Pr 



(2) 



By differentiating both expressions and solving for the 
time when the signals are minimized, one may write: 



+ T = 



-7 + ipivot 

1 



(3) 



Differentiating both RV models with respect to time 
twice and evaluating at the moment when both signals 
are minimized, yields: 



Kc 
P 2 



7 
4tt 2 



(4) 



Equations may therefore be used to determine 
some information about HAT-P-31c. 

To evaluate the statistical significance of HAT-P-31c, 
we performed an F-test between the one-planet and two- 
planet models. In both cases, HAT-P-31b is assumed to 
maintain non-zero orbital eccentricity. Assuming HAT- 
P-31c is on a circular orbit, the false-alarm-probability 
(FAP) from an F-test is 3.0 x 10 -12 , or 7.0-cr. Assuming 
HAT-P-31c is on an eccentric orbit requires 2 more de- 
grees of freedom and thus reduces the FAP to 1.3 x 10~ 10 , 
or 6.4-cr. From a statistical perspective then, the pres- 
ence of HAT-P-31c is highly secure. We point out this 
determination of course assumes purely Gaussian uncer- 
tainties and no outlier measurements. 

3.2.3. Fitting algorithm 

We utilize a Metropolis-Hastings Markov Chain Monte 
Carlo (MCMC) algorithm to globally fit the data, in- 
cl uding the stellar density p rior (our routine is described 
in iKipping fc Bakos I (|2011[ )). To ensure the parameter 
space is fully explored, we used 5 independent MCMC 
fits which stop once 1.25 x 10 5 trials have been accepted 
and burn-out the first 20%. This leaves us with a total of 
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TABLE 6 

Comparison of RV models attempted for HAT-P-31 

Model x 2 BIC a 

Circular Planet 3714.9 3729.6 

Eccentric Planet 204.2 226.3 

Eccentric Planet + Drift 159.6 185.5 

Eccentric Planet + Trojan 191.1 216.9 

Eccentric Planet + Drift + Trojan. . . 113.2 142.7 

Eccentric Planet + Drift + Quadratic 33.9 63.4 

Eccentric Planet + Circular Planet . . 33.6 66.8 

Eccentric Planet + Eccentric Planet . 33.2 73.8 

a BIC — Bayesian Information Criterion (Schwarz 1978; 
ILiddlell2007l) 

5 x 10 5 points for the posterior distributions. At the end 
of the fit, a more aggressive downhill simplex x 2 mini- 
mization is used, for which the final solution is used for 
Figures 0^3] 

There were 14 free parameters in the global fit: {r b , p 2 , 

[Ti] b ,b b ,P b , e b sinw b , e b cosw b , 7 Kc ck, 7fies, 7Subaru, Kb, 
7, 7, OOT}, which we elaborate on here. T b is the time 
of transit minimum (Kipp mg"ll2011l) . frequently dubbed 
by the misnomer "mid-transit time", p 2 is the ratio-of- 

radii squared and P b is the orbital period. \Ti\ b is the 
"one-term" approximate equation (Kip ping! I2010D for 
the transit duration between the instant when the center 
of the planet crosses the stellar limb to exiting under the 
same condition. b b is the impact parameter, defined as 
the sky-projected planet-star separation in units of the 
stellar radius at the instant of inferior conjunction. e b 
is the orbital eccentricity and uj b is the associated po- 
sition of pericenter. 7 terms relate to the instrumental 
offsets for the radial velocities. Similarly, OOT is the 
out-of-transit flux for the HATNet photometry. Finally, 
K b is the radial velocity semi-amplitude and 7 (drift) & 
7 (curl) are the first and second time derivatives of 7. 

Final quoted results are the median of the marginalized 
posterior for each fitted parameter with 34.15% quan- 
tiles either side for the 1-cr uncertainties (see Table [7]). 
The uncertainties on p, p 2 and Rp have been conserva- 
tively doubled for reasons described in [J2] Histograms of 
the posterior distributions for the fitted parameters are 
provided in Figure [5j We find that the stellar jitter of 
this star is at or below the measurement uncertainties 
(~ 2m/s). 

Table [7] provides estimates for some minimum limits 
on various parameters of interest relating to HAT-P-31c. 
These limits are determined by the known constraints on 
the minimum P c . We determined this value by forcing 
a circular orbit Keplerian fit for planet c through the 
data, stepping through a range of periods from 1 year 
to 5 years in 1 day steps. The minimum limit on P c is 
defined as when A^ 2 = 1, relative to the quadratic trend 
fit, occurring at 2.8 years. 

4. DISCUSSION 

4.1. Physical Properties of HAT-P-31b&c 

HAT-P-31b is a 2.17 Mj hot- Jupiter transiting the host 
star once every 5.005 days. Due to the lack of follow-up 
photometry obtained for this object (as a consequence of 
the nearly integer orbital period) , we have only HATNet 
photometry, which is of lower signal-to-noise than dedi- 
cated follow-up. This fact, combined with our choice to 



TABLE 7 
Global fit results for HAT-P-31 

Parameter 3 Value 
Fitted Parameters 

n [days] 5.005425^™^ 

r b [BJDtdb - 2,450,000] 4320.8866±o ^ 

[fi) b [ S ] 16300fS» 

p 2 b l%\ o.65±g;li 

h 0.57i°;g 

oot i -nnotii •;;;!;!;!;!;; 

K b [ms" 1 ] 232.51" 

e b sinu b -0.2442±°;°°g 

e b coso; b 0.0185^;°°?° 

7Kcck [ms" 1 ] 29.0 ' ! ; j 

TSubaru [ms" 1 ] 18-8l|;l 

7FIES [ms" 1 ] 92.711° 

7 [ms-May- 1 ] 0.141±g;gg 

7 [ms-^ay 2 ] 0.002261°;°°°^ 

SME Derived Quantities 

ui be 0.2078* 

u 2 hc 0.3550* 

p* CB [gem- 3 ] 0-69tg-Je 

HAT-P-31b Derived Properties 

w i> (K)(14 

pt. 2450+ - 0045 

Ub[°] 274.3ti'j 

log( Sb [cgs]) 3.61±°;^ 

Mr*) 8.9 • \i 

ibl°] 87.lj£? 

[TiAh W 18500™ 

[T2,s] b w I'^oo 

p» ".080 

Mb [Mj] 2.171±g"J?? 

R b [Rj] • ov;:;v,: 

Corr(R b ,M b ) 0.795 

Pb [gem- 3 ] 2.18±J || 

a b [au] iui55 •;; !;!:; 

[T cq ] b [K] 1450l 2 ?° 

<■->■ 0.19Qt ;°i 

F pcri [lO^rgs-icm" 2 ] 1.69±J |2 

F ap [lOSergs-icm- 2 ] 0.62±°- % 

(F) 1 [lO^rgs^cm- 2 ] 0.99tgj 2 

HAT-P-31c Derived Quantities 

t c + 0.25P C [BJD TDB - 2,450,000] 5254.9+g ! 

K C P~ 2 [ms-May- 2 ] 7.64±°-^ 

P c [years] > 2.8 

K c [ms" 1 ] > 60 

M c [Mj] > 3.4 

a Quoted values are medians of MCMC trials with errors 

given by 1-ct quantiles. "b" subscripts refer to planet HAT- 

P-31b and "c" subscripts refer to HAT-P-31c. 

k Fixed parameter 

c Parameter is treated as a prior 

^ Parameter determined using SME and YY isochrones 

e The Safronov number is given by — ■^(V csc/Kprb) 2 — 

(a/R p )(M p /M ir ) (see[H ansen Barman||2007D . 

^ Incoming flux per unit surface area, averaged over the 

orbit. 
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Fig. 5. — Posterior distributions of the fitted parameters used in the global fits (described in i]3.2[l from our global fits. Histograms 
computed from 5 X fO 5 MCMC trials. Bottom-right panel shows joint-posterior of the radial velocity drift and curl, with white denoting 
the 2-ct region of confidence and gray the 3-cr. 



double all uncertainties on depth related transit terms, 
leads to a large uncertainty on the planetary radius of 
Rb = 1.07_to*32 R-Ji consistent with many other hot- 
Jupiter objects (see http://exoplanet.eu). 

High-precision radial velocities also indicate the pres- 
ence of an outer body, HAT-P-31c, found through an 
induced quadratic trend in the RV residuals. Keplerian 
fits are unable to convincingly distinguish between a cir- 
cular or eccentric orbit for this object. HAT-P-31c has 
a minimum mass of M c > 3.4 Mj and eccentric orbit 
solutions significantly increase this figure. It is unclear 
whether HAT-P-31c is a brown dwarf or a "planet" , and 
future work will be needed to determine this. 



4.2. Orbital Stability 
4.2.1. Circular fit for HAT-P-Slc 

Here we discuss our procedure to test the dynamical 
stability of two possible orbital configurations. It should 
be noted that the eccentricity of HAT-P-31b is highly 
secure but the eccentricity of HAT-P-31c remains un- 
clear. For this reason, we repeat our simulations assum- 
ing both a circular and eccentric orbit for HAT-P-31c, 
begi nning with the former. We utilize the Systemic Con- 
sole (jMeschiari et al. Il2009f) for this purpose assuming a 
coplanar configuration. Employing the Gragg-Bulirsch- 
Stoer integrator, orbital evolution was computed for 
250,000 years for the HAT-P-31 system (see Figure [7]). 
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Fig. 6. — Top row: Marginalized posterior distributions for K c and P c when we attempted to fit for a second Kcplcrian signal, instead 
of a quadratic trend. The multi-modal nature of these histograms reflect the unconverged nature of the fits. Bottom row: History of the 
MCMC trials for the same parameters and the same fit. Each continuous line represents one of the ten independent MCMC chains. The 
lines clearly illustrate the inability of the current data to converge upon a solution for HAT-P-31c. 



We first consider the circular case. The orbital period 
and mass of HAT-P-31c are non-convergent parameters 
and so we can only provide an orbital solution which 
gives a good fit to the data, but is not necessarily unique. 
To this end, we proceeded to input HAT-P-31c with 
P c = 4.86 years and M c — 13.0 Mj, corresponding to the 
solution presented in Table [51 This test revealed minor 
evolution over the course of our simulations, indicating 
a stable and essentially static configuration. 

4.2.2. Eccentric fit for HAT-P-31c 

To test the eccentric fit, we again used the lowest 
X 2 solution presented earlier in Table [SJ corresponding 
to P c = 4.82 years and e c = 0.285. We found that 
the system was also stable over 250,000 years (see Fig- 
ure [7]). However, the simulations do show the eccentric- 
ity of planet b varying sinusoidally over a timescale of 
~125,000years with an amplitude of ~ 0.01. The eccen- 
tricity of planet c also varies in anti-phase but with a 
much smaller amplitude. The eccentric evolution shows 
faster apsidal precession for planet b, but this is unlikely 
to be observable through changes in the transit dura- 
tion. We estimate the duration will change by 0.2 s over 
10 years (co rresponding to A cu b = 0.027°) using the ex- 
pressions of IKipping I pOloh . 

The orbital period and semi-major axes of both bodies 
were stable over the 250,000 years of integration consid- 
ered here. 

4.2.3. Habitable-zone bodies 

We tried adding a habitable-zone Earth-mass planet 
on a circular orbit into the system and testing stabil- 



ity. One may argue that the probable history of this 
system involved the inwards migration of HAT-P-31b 
and that this migration through the inner protoplane- 
tary disk would essentially eliminate the possibility of an 
Earth -like planet forming i n the habitable-zone. How- 
ever, iFogg fc Nelson! ([2007) have shown that this not 
necessarily true. In their simulations, it is found that 
> 60% of the solid disk survives, including planetesimals 
and protoplanets, by being scattered by the giant planet 
into external orbits where dynamical friction is strong 
and terrestrial planet formation is able to resume. In 
one simulation, a planet of 2 M s formed in the habitable- 
zone after a hot- Jupiter passed through and its orbit sta- 
bilized at 0.1 AU. 

For a planet to receive the same insolation as the Earth, 
we estimate P — 604 days. For our circular orbit solution 
of HAT-P-31c, the habitable-zone Earth-mass planet is 
stable for over 100,000 years. For our eccentric orbit so- 
lution, the Earth-like planet is summarily ejected in less 
than 1000 years. 

4.3. Circularization Timescale 

Due to the poor constraints on the planetary radius, 
there is a great deal of uncertainty in the circulariza- 
tion timescale (T c j rc ) for HAT-P-31b. Neverth eless, us- 
ing the equations of I Adams fc Laughlinl (|2006l ). we used 
the MCMC results to compute the posterior distribution 
of rare- We find that the age of HAT-P-31 is equal to 
24iis 2 circularization timescales, assuming Qp = 10 5 . 
This indicates that we currently have insufficient data 
to assess whether the observed eccentricity is anomalous 
or not. Improved radius constraints will certainly aid in 
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Fig. 7. — Two possible realizations for the orbital evolution of planets HAT-P-31b and c. The solid lines show the evolution starting 
from an eccentric orbit solution for HAT-P-31c. The dashed lines show the evolution starting from a circular orbit solution for HAT-P-31c. 



this calculation and may lend or detract credence to the 
hypothesis of eccentricity pumping of the inner planet by 
HAT-P-31c. 
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